This invention is relevant to seismic data processing in the field of geophysical exploration for petroleum and minerals. The general seismic prospecting method involves the transmission of elastic, or "seismic," waves into the earth and reception of reflected and/or refracted waves at the earth's surface (or, occasionally, in a wellbore) via geophones, hydrophones, or other similar devices (hereinafter referred to collectively as "geophones"). The elastic waves are generated by sources such as dynamite, air guns, and hydraulic vibrators. As these waves propagate downward through the earth, portions of their energy are sent back to the earth's surface by the acts of reflection and refraction which occur whenever abrupt changes in impedance are encountered. Since these impedance changes often coincide with sedimentary layer boundaries, it is possible to image the layers by appropriate processing of the signals returned to geophone groups.
Raw seismic data as they are recorded are not readily interpretable. While they show existence of formation interfaces, raw data do not accurately inform the interpreter as to the location of these interfaces. The process of migration, also called imaging, is well known in the art, and is the repositioning of seismic data so that a more accurate picture of subsurface reflectors is given. In order to perform migration calculations, the seismic velocities of the subsurface at a multitude of points must first be determined or approximated. In the discussions to follow, the unqualified term "velocity" will be used as a short-hand expression that means the velocity of propagation of an acoustic wavefield through earth formations. Generally, for the purposes of this discussion, the term is meant to apply to the propagation of compressional waves although shear waves are not excluded. The term "velocity model" is used to refer to the two- or three-dimensional spatial distribution of velocity. The large-scale velocity model covering the extent of the seismic data acquisition volume is usually a complicated structure with vertically and laterally varying velocity. The commonly used Kirchhoff migration method is well known in the art and requires the calculation of traveltimes through the velocity model. "Traveltime" means generally the amount of time a seismic signal takes to travel from seismic source to a subsurface interface reflection point to a seismic receiver. Current methods of computing traveltimes necessary for two-dimensional and three-dimensional depth migration and associated velocity estimation are inefficient or potentially error-prone when applied to the complicated large-scale velocity models typically encountered.
There are currently two general methods of computing the grid of traveltimes needed to migrate data: the ray tracing methods, and the methods which seek the direct solution of the eikonal equation. As is known to one skilled in the art, the eikonal equation is a form of the wave equation for harmonic waves, valid only where the variation of properties is small within a wavelength, otherwise termed "the high-frequency condition." The eikonal equation relates the gradient of traveltime to the velocity structure.
Seismic ray tracing methods are implemented by solving the ordinary linear differential equations which are derived by applying the method of characteristics to the eikonal equation, a technique known to those skilled in the art. Ray tracing allows the determination of arrival times throughout the subsurface, by following raypaths from a source location, which raypaths obey Snell's law throughout the subsurface volume for which the velocity model is known. Traveltimes along the rays are then interpolated onto a grid of the subsurface. The ray equations may be solved with "shooting" methods or with "bending" methods, both of which are well known to those skilled in the art.
Shooting formulates the ray tracing equations into an initial-value problem, where all ray direction and position components are defined at the source location at time zero. Then equations are recursively solved to trace the rays throughout the medium. Bending formulates the ray tracing equations into a two-point boundary value problem and is based on Fermat's principle, which states that the seismic raypath between two points is that for which the first-order variation of traveltime with respect to all neighboring paths is zero, and attempts to locate a raypath between two points by determining a stationary traveltime path between them. Shooting is generally more efficient computationally than bending; however, both approaches present difficulties and potential inaccuracies when used to compute the gridded traveltime fields required by two- and three-dimensional depth migration. These difficulties and inaccuracies increase with the degree of complexity of the velocity model and with the increase in overall size of the large-scale velocity model. Depth migration typically requires robust grids of traveltimes for high quality images. As used herein, the term "robust" means a process which reliably generates accurate grids of traveltimes regardless of the complexity of the large-scale velocity model.
Complications arise in shooting and bending calculations, especially in large-scale complex velocity models, because each ray is computed independently of all others, and because small changes in the angle of incidence may lead to large changes in ray direction. The latter complication is a manifestation of the nonlinearity of the problem of ray tracing. Thus, with shooting methods, it is difficult to obtain the uniform ray coverage throughout the large-scale model that will admit traveltime interpolation onto a grid. Also, rays may cross, indicating the natural multivaluedness of wave propagation. The multivaluedness of the wave propagation increases in proportion to the complexity and overall size of the velocity model. The multivaluedness of the wavefield complicates the computation because it is desirable to track the multivalued wavefronts, and it complicates the gridding because one may not interpolate between different branches of a multivalued wavefront.
Ray shooting itself is quite efficient, but a robust application using it to compute grids of traveltimes may not be efficient for large-scale velocity models exhibiting complexities which typically give rise to multivalued traveltimes. While the bending method will determine a raypath to any point in the grid, it is subject to local, instead of global, minimization, resulting in traveltime errors in large-scale models. Bending is also inefficient for determining three-dimensional grids of traveltimes.
In large-scale velocity models, crossing ray trajectories and refractions cause complicated amplitude and phase behavior of seismic wavefields. In order to accurately image seismic data, it is necessary to account for these amplitude and phase effects. The calculation of these amplitudes and phase changes along ray trajectories adds significantly to the computational cost of the methods.
Because of their efficiency, methods based on direct solutions of the eikonal equation are particularly attractive for the calculation of migration traveltimes. The calculation is performed by finite-differencing the eikonal equation on either a cartesian grid, or on a polar/spherical grid which is readily interpolated to a Cartesian grid. The traveltimes are piecewise smooth and they fill the computational grid. While methods based on finite-differencing the eikonal equation are efficient and computationally well suited for migration, the solution calculates so called first-arrival traveltimes.
In a large-scale depth model where there are considerable relative changes in velocity, these first-arrival traveltimes often correspond to non-energetic portions of the seismic wavefield. This creates problems for migration since it is critical that the traveltimes used in migration correspond to the energetic portions of the wavefield. Changes in velocity distort the shape of the wave propagation front and create opportunities for frequency components to separate, for headwaves to develop, and for triplications to occur. When high velocity zones are present, the first-arrivals may be non-energetic headwaves. As energy propagates in complicated large-scale models, raypaths tend to eventually cross, which causes phase shifts and triplications. First-arrival traveltimes follow the fastest branch of the triplication bow-tie, which is also the low energy branch. Traveltimes calculated in a large-scale velocity model with a finite-difference method are susceptible to these types of inaccuracies and cannot accurately parameterize the asymptotic Green's functions required for Kirchhoff imaging.
Kirchhoff migration is generally accepted to be the most efficient and practical method of imaging two- and three-dimensional prestack seismic data. However, Kirchhoff algorithms using first-arrival traveltimes typically do a poor job of imaging complex structures, as observed for example by Gray and May, Geophysics 59: 810-817 (1993), and Geoltrain and Brac, Geophysics 58: 564-575 (1993). Even ray tracing methods which calculate multiple arrivals and most energetic arrivals along with estimates of amplitude and phase do not always result in satisfactory images, as indicated by Audebert et. al., (1995) "Imaging Complex Geologic Structure with Single-Arrival Kirchhoff Prestack Depth Migration," scheduled for publication in Geophysics 62, No. 5 (1997).
It is generally accepted that migration methods which use recursive wavefield continuation to backwards propagate the received wavefield produce the best images. These migration methods, commonly referred to as "finite-difference migration" or "phase-shift migration" are well known to those skilled in the art. The problems with recursive finite-difference and phase-shift methods are that they require regular spatial sampling, are computationally intensive, and are not readily applicable to three-dimensional prestack migration. This makes the methods practically applicable only to two-dimensional regularly sampled data, severely limiting their utility. Methods based on the Kirchhoff integral are thus attractive, especially for three-dimensional prestack imaging objectives. Kirchhoff algorithms can easily accommodate irregular sampling and they can be applied in a target-oriented fashion. Using first-arrival traveltimes is popular because they are efficiently computed and they have the attractive property of filling the entire computational grid.
In their 1993 Geophysics article, Geoltrain and Brac ask the question "Can we image complex structures with first-arrival traveltime?" They conclude that they encounter imaging difficulties with Kirchhoff migration using first-arrival traveltimes, and propose to either ray trace to find the most energetic arrivals, or to calculate dynamically correct multiple-arrival Green's functions. Nichols ("Imaging Complex Structures Using Band Limited Green's Functions," Ph.D. Thesis, Stanford University, 1994) calculates band-limited Green's functions to estimate the most energetic arrivals. He estimates not only traveltime, but also amplitude and phase. Both these traveltime calculation methods are computationally complex and much more costly than first-arrival traveltime computation methods.